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Ring-like ripples on the surface of icicles are an example of morphological instability of the ice- 
water interface during ice growth under supercooled water film flow. The surface of icicles is typically 
covered with ripples of about 1 cm in wavelength, and the wavelength appears to be almost inde- 
pendent of external temperature, icicle radius, and volumetric water flow rate. One side of the 
water layer consists of the water-air surface and growing ice is the other. This is one of the more 
. complicated moving phase boundary problems with two interfaces. A recent theoretical work [K. 

Ueno, Phys. Rev. E 68, 021603 (2003)] to address the underlying instability that produces ripples 
is based on the assumption of the absence of airflow around icicles. In this paper, we extend the 
previous theoretical framework to include a natural convection airflow ahead of the water-air surface 
and consider whether the effect of natural convection airflow on the wavelength of ripples produced 
on an ice surface is essential or not. 
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I. INTRODUCTION 



Little is known on the study of morphological instability of the solid-liquid interface when a thin layer of moving 
fluid separates the developing solid from its surrounding. For example, the problem of icicle growth involves complex 
moving boundary problems with phase change. When an icicle grows, a thin water film from the melting snow and 
ice at the root of the icicle flows down along its surface and refreezes onto it by releasing latent heat of solidification 
to the ambient air below °C. During the icicle growth, ice does not grow uniformly, but ring-like ripples are often 
observed on its surface. [l[ By supplying water continuously from the top of a wooden round stick and of a gutter 
on an inclined plane set in a cold room below °C, a ripple pattern similar to that observed on natural icicles is 
produced on the ice surface. Surprisingly, the distance between two peaks of ripples experimentally produced as 
well as that of natural icicles always measures around a centimeter scale. 

Theoretical works aimed at explaining the underlying dynamic instability that produces ripples are recent. [Mi a 
stability analysis for the ice- water interface disturbance was developed based on heat flow in the water and atmosphere, 
and thin film water flow dynamics. From the initial model, it was found that the ripple wavelength is determined 
from A = 27r/io-Pe;/a ma x, and that the ripples should move down the icicle. [3j Here ho is the mean thickness of the 
| water layer, Pei is the Peclet number, which is a dimensionless number defined as the ratio of the heat transfer due 
to the water flow to that due to the thermal diffusion in the water layer, and a max is a dimensionless wave number 
at which the amplification rate of the ice-water interface disturbance acquires a maximum value. By considering 
different boundary conditions from those used in the initial model, a quite different ripple formation mechanism was 
developed. [1, Q A new formula to determine the wavelength of ripples was derived: A = 2Tr(a 2 hoPei/3) 1 ^ 3 , which 
contains two characteristic lengths ho and a. Q Here a is the capillary length associated with the surface tension of 
the water-air surface. In the new model, the influence of the shape of the water-air surface on the growth condition of 
the ice-water interface was taken into account. Therefore, another length scale a was introduced. The new model also 
predicted that ripples should move upward. The upward ripple translation was already suggested by the observation 
that many tiny air bubbles were trapped in the upper side of any protruded part of ripples during the icicle growth, 
and lined up upward. [l| However, there was no theoretical explanation for the upward ripple translation mechanism. 

Both models yield one-centimeter scale wavelength, but the translational direction of the ice ripples is opposite. 
Recently we solved numerically the same governing equations with the same boundary conditions as those used in the 
initial model. However, the numerically obtained amplification rate of the ice-water interface disturbances showed 
positive values for all wave numbers, Q which means that a max does not exist and there is no mechanism to select a 
characteristic length. On the other hand, the analytical results for the amplification rate and the translation velocity 
of ice ripples obtained in the new model were in good agreement with those numerically calculated. Moreover, there 
was also good agreement between the theoretical predictions of the dependence of ripple wavelength on slope angles of 
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FIG. 1: Schematic view of the layer ahead of the water-air surface. Ice is covered with a supercooled water film. The x axis is 
parallel to the direction of the supercooled water flow and the y axis is normal to it. Ti a is the temperature at the water-air 
surface, (a) is the situation in absence of airflow. A linear air temperature distribution T a (y) was assumed, (b) is the situation 
in presence of airflow. U a (x, y) and T a (x, y) are undisturbed velocity and temperature distributions, ho and 5 are the thickness 
of the water layer and that of the thermal boundary layer, respectively, g is the gravitational acceleration and 9 is the angle 
with respect to the horizontal. The flowing supercooled water layer, not to scale, is much thinner than the thickness of the 
thermal boundary layer. 



the inclined plane and water supply rates and our experimental results. Finally, upward ripple motion at about half- 
speed of the mean growth rate of icicle radius was observed experimentally as theoretically predicted, but downward 
traveling ripples were not observed. [7| 

In the previous theoretical models, [J-Hl ice was covered with a supercooled water layer and there was no airflow 
around icicles. The latent heat released at the ice-water interface was assumed to be transferred in the air by thermal 
diffusion through the water layer. As shown in Fig. [T](a), for simplicity, a linear air temperature distribution T a (y) 
was assumed. j5| From the energy conservation at the ice- water interface and water-air surface, the mean growth rate 
of icicle radius is given by V — —K a Too/(LS), where K a is the thermal conductivity of air, L is the latent heat per 
unit volume and Too is the air temperature at a distance S from the water-air surface. || The water layer of thickness 
ho changes by varying the water supply rate. [sj. [Io| However, since V does not include ho, the icicle growth rate 
does not depend on the water supply rate. Since V contains the parameters Too and 5, the ice growth rate 

is controlled by the rate of latent heat loss from the water-air surface to the surrounding air. However, it was not 
possible to estimate the value of V because the physical meaning of the assumed distance 6 was unclear. 

Recently, the growth of icicles has been treated as a free boundary problem to find an ideal growing shape for 
icicles. Q The latent heat transferred from the icicle surface to the surrounding air through the water layer leads 
to an increase in air temperature and to a change in density because it is temperature dependent. If the density 
decreases with increasing temperature, buoyancy force arises, and warmer air moves up along the ice surface. This 
effect is restricted to a thin layer ahead of the water-air surface, as shown in Fig. [T] (b). Short et al. emphasized the 
importance of heat transfer through such a convective boundary layer around icicles, and derived a formula for the ice 
growth velocity normal to the icicle's surface. The form is the same as V mentioned above, but a critical difference 
is that the length 5 in paper [8] is the boundary layer thickness. Hence, it was possible to estimate the values of S 
and V if the value of an unknown parameter in S was given. 0, Q It was also suggested that similarity solutions 
for the coupled Navier-Stokes and heat transfer equations in the Boussinesq approximation can provide the basis for 
understanding of the boundary layer. [f| In this paper, the value of the unknown parameter in 5 is determined by 
obtaining the similarity solutions. 
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Since heat transfer can be greatly influenced by the upward natural convection airflow, the question is whether the 
enhancement of heat transfer due to convection affects the wavelength of ripples on icicles. On the other hand, it 
is known that the wavelengths are almost independent of the length of the icicles and the ambient air temperature. 
In order to clarify these problems, in this paper, a linear stability analysis was performed on the ice-water interface 
disturbance during the ice growth in the presence of a supercooled water film flow and a natural convection airflow. 

II. THEORY 

Instead of dealing with the elongated carrot-shaped geometry of the icicle, H ice growth on a flat gutter on an 
inclined plane of finite length will be considered. The following theoretical analysis is restricted to two-dimensional 
vertical cross-sections of the gutter, as shown in Fig. [1] The origin of the x axis is the bottom of the gutter and the y 
axis is normal to it. What is new here is that the effect of a natural convection airflow is being incorporated into the 
previous theoretical frameworks 0-0] with modifications of some of boundary conditions, letting us treat synthetically 
heat flow in the ice, water and air through a disturbed ice- water interface and water-air surface, as well as thin water 
film flow and airflow. 

A. Governing equations 

The velocity components in the x and y directions in the water layer, ui and vi, are governed by the Navier-Stokes 
equations driven by gravity and the continuity equation: 10] 

dui dui dui 1 dpi ( d 2 ui d 2 ui \ 

' u i^z + v i^r: = -- -sz+^i uzt + itt -gsmB, (l) 



dt dx dy pi dx \ dx 2 dy 

dvi | ^dvi | ^dyi _ 1 dpi | ^ / d 2 v t | d 2 vi \ gCQZ Q i 2 ) 
dt dx dy pi dy \ dx 2 dy 2 J 

dui ^ dvi q 
dx dy 

where v% = 1.8 x 10~ 6 m 2 /s and pi = 1.0 x 10 3 kg/m 3 are the kinematic viscosity and the density of water, g the 
gravitational acceleration, pi the pressure in water. 9 is the angle with respect to the horizontal, as shown in Fig. [T] 
On the other hand, employing the Boussinesq approximation, the velocity components in the x and y directions in 
the air, u a and v a , are governed by the following equations driven by buoyancy force and the continuity equation: [Toj 

du a du a du a 1 d(p a -p a0 ) ( d 2 u a d 2 u a \ . 

+ u a^~ + v a — = — + v a -— + -— + gp{T a - Too) sin 5, (4) 



dt dx dy poo dx \ dx z dy 

9v a dv a dv a 1 d{p a -p a0 ) (d 2 v a d 2 v a \ 

-zrr + u a-^~ +v a -— = 5 h v a -7T^- + -z-T +9P{ T a - Toojcosfl, (5) 

dt dx dy poo dy \ dx z dy z J 

dUa , dVa n ra\ 
^ + ^ = °' (6) 

where p a is the pressure in air, p a Q the static pressure, poo the density of air at the temperature Too, v a = 1.3 x 10~ 5 
m 2 /s and /3 — 3.7 x 10~ 3 K^ 1 are, respectively, the kinematic viscosity and the volumetric coefficient of thermal 
expansion for air. The continuity equations © and ([5]) can be satisfied by introducing the stream functions ipi and 
ipa such that ui = dipi/dy, vi = —dipi/dx, u a = dip a /dy and v a = -dip a /dx. 

Neglecting viscous dissipation in the energy equation, the equations for the temperatures in the ice T s , water T; 
and air T a are [To| 

3T, + (7) 



dt \ dx 2 dy 2 
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dTi | &Ti | v m _ r fd 2 Ti | d 2 ^ 
<9t c?a: <9?/ \ <9x 2 <9z/ 2 



ar a , dr a , 9r a ^a 2 T a , 9 2 r a 



ot ox dy \ ox z dy 



(9) 



where k s = 1.15 x 1(T 6 m 2 /s, m = 1.33 x 10 -7 m 2 /s and n a = 1.87 x 1CT 5 m 2 /s are the thermal diffusivities of 
ice, water and air, respectively. Equations (@|, ©, @ and © are new part that has been added to the previous 
formulation. 0-0] 



B. Boundary conditions at the ice-water interface and water-air surface 

1. Hydrodynamic boundary conditions 

Neglecting the density difference between ice and water, both velocity components ui and vi at a disturbed ice- water 
interface, y = £(t, x), must satisfy the no-slip condition: [3j 



ui\ y=c = 0, vi\ y=< =0. 
The kinematic condition at a disturbed water-air surface, y = £(t,x), is @ 

at dt 

The continuity of velocities of water film flow and airflow at the water-air surface is 

Ul\y^ = U a \ y = i , Vl\ y=i = V a \ y=i . 

The condition for continuity of shear stress at the water-air surface is [HI 

Pin 



(10) 



(11) 



(12) 



(13) 



The difference of the normal stress on either side of the water-air surface must be the capillary force resisting dis- 
placement: [Iol - fl2| 



/ dui 




dvi 




( du a 


dv a 




V dy 


H 

y=i 


dx 






v=t dx 


y=t) 



I l o dVa 

■ Pa\y=(, + ZPaVa^r- 

dy 



dvi 

-Pi\ y=( + 2pivi — 



eft 
' dx 2 



(14) 



where 7 = 7.6 x 10~ 2 N/m is the surface tension of the water-air surface. The boundary conditions (fTCTj) and (jTTj) 
are the same as those used in the previous papers. 0-0| Since an airflow is taken into account in this paper, the 
continuity condition of the water film and airflow velocities at the water-air surface is a new part, and the shear and 
normal stress conditions are modified from those in the previous papers. 043 



2. Thermodynamic boundary conditions 



The following boundary conditions are exactly the same as those in the previous papers. 043] The continuity 
condition of temperature is imposed at the ice-water interface: 



T i\y=C - T t 



s\y=C 



T s i + AT. 



si •■ 



(15) 



where T s i is the temperature at the flat ice-water interface and AT S ; is a deviation from it when the ice- water interface 
is disturbed. The energy conservation at the ice-water interface is 



K, 



dT\ 
dy 



y=C 



(16) 
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where L = 3.3 x 10 8 J/m 3 is the latent heat per unit volume, and K s = 2.22 J/(mKs) and Ki = 0.56 J/(mKs) are 
thermal conductivities of ice and water, respectively. 

The continuity condition of temperature is imposed at the water-air surface: 

T l\v=£ = T a\y=£ = Tia, (17) 

where 7] a is a temperature at the water-air surface. The energy conservation at the water-air surface is 

dTi 



Kl 7T 
ay 



v=t ay 



(18) 



where K a — 0.024 J/(mKs) is the thermal conductivity of air. 

In the initial model, Q the continuity condition of temperature at the ice- water interface and water-air surface was 
T s \ y= £ = Ti\ y= Q — T s i and Ti\ y= ^ — T a \ y= ^ = T; a + A7] a , where A7] a is a deviation from Tj a when the water-air 
surface is disturbed. Instead, we use Eqs. (ITBI) and (|17l) as in the previous papers. [3-0] The difference in these 
boundary conditions led to critically different results between two models mentioned in the Introduction in this 
paper. When the chemical potential of water equals that of ice, it seems reasonable to assume the boundary condition 
T s \ y= £ = Ti\ v= q = T s i at the ice- water interface, then T s i is the equilibrium freezing temperature (T s i = °C for 
pure water). As in paper, [l]| however, the chemical potential of water is not necessarily equal to that of ice because 
the ice-water coexistence considered here is expected to be in a non-equilibrium state in the presence of external 
disturbance at the water-air surface and shearing water flow. The deviation AT s i at the ice- water interface caused by 
the external disturbance does not disappear by thermal diffusion in the water because the thermal relaxation time 
for the temperature fluctuation with about 1 cm corresponding to ice ripple wavelength is much longer than the time 
defined by the inverse of the shear rate of water film flow considered here. In other words, the equilibrium state at 
the ice- water interface is not attained in the presence of shearing water flow. [7j We will see in IIIIDI that AT 5 j is 
dependent on the temperature distribution in the water layer subject to the external disturbance at the water-air 
surface. On the other hand, since shear stress has a value of zero at the water-air surface, the deviation AT; Q at the 
water-air surface disappears by thermal diffusion in the air. Hence, the temperature at the water-air surface remains 
at 7} a , which will be determined in IIII Al 



C. Perturbation 



As shown in Fig. [TJ only a one-dimensional perturbation in the x direction of the ice-water interface with a small 
amplitude £fc is considered: C,(t,x) = (fcexpfcrf + ikx], where k is the wave number and a = + icrW. Here <j( r ) 
and v p = — o-W/fe are the amplification rate and the phase velocity of the perturbation, respectively. ipi^ ipa-> Pi-, 
p a , T s , T\ and T a are separated into unperturbed steady and perturbed parts as follows: £ = ho + tfti = ipi + tpL 
V> = + V4, PI = Pl+ P'i, Pa = P a + p' a , T s = T s + 7* T ; = T ( + T{ and T a = T a + T' a . The corresponding 
perturbation of the water-air surface with a small amplitude is t;'(t,x) = £fcexp[ert + ikx]. As in the previous 
papers, 0-0] the following calculation is based on a linear stability analysis taking into account only the first order 
of £fc . The quasi-stationary approximation is also used: the time dependence of the perturbed part of equations can 
be neglected because the time evolution of the ice- water interface perturbation is considerably slow compared to that 
of the above perturbation fields. [l4j 



D. Equations of flow and temperature distributions in the air boundary layer 

Natural convection airflow considered here are restricted to a boundary layer regime and to conditions that lead to 
a similarity solution, that is, to a description of the flow by ordinary differential equations and boundary conditions 
in terms of a single coordinate r](x,y). Under this assumption the unperturbed quantities ij) a {x,y) and T a (x,y) are 
expressed as follows :[l5| 

- - T — Too 

ipa = UaoS F a (ri) = v a GrF a (rf), T a * = a - °° , (19) 

-L la 1 oo 

where 77 = (y — h )/5 , S — 4x/Gr and u a0 = v a Gr 2 /(4x). Here Gr = 4(Gr 2 ,/4) 1 / 4 is the modified local Grashof 
number, Gr x — g(3AT a x 3 being the local Grashof number. AT a = T\ a — T x is the temperature difference between 
the water-air surface and the ambient air temperature far away, x is the distance from the bottom of the gutter. 
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Applying the boundary layer approximation to the Boussinesq equations (j4]), (|5j), {[6|) and ©, t^ a (x, y) and T a (x, y) 
are governed by [TH [l6j 

d$a d 2 ^ a dtp a d 2 ^ a d 3 ^ a - . 

-5 — jr~^ « o-o- =^0-5-5- +gp{T a -T 00 )sm6, (20 

ay oxoy ox ay 2 - ay 6 

C^adTa _ Chkdfa = ^ <Pfa 

dy dx dx dy a dy 2 

When Eq. (|TT))) is substituted into Eqs. (j2T))) and (|2"Tj) . the dimensionless functions F a and T a * are obtained from the 
two coupled ordinary differential equations: [1 0| 

^ = - 3f # +2 ( f |) 2 - f " si " e ' (22) 
" - 3Pr - p °% (23 » 

where Pr a = v a / n a — 0.7 is the Prandtl number of air. 

We assume stream function disturbance ip' a and temperature disturbance T' a in the air to be of the form: 

V4 = u a0 fa(v)^kexp[at + ikx), T' a = H a (-q)G a £, k exp[at + ikx), (24) 

where f a and H a are the dimensionless disturbance amplitude functions, and G a = —dT a /dy\ y —h . When ip a = ip a +ip' a 
and T a — T a + T' a are substituted into the complete equations Q , ([5]) and © , we obtain the differential equations 
for the functions f a and H a : 



d\fa 

drf 



= -3F Q 



drf 



dF a \ d 2 f a 



drj J drj 2 



14 3F a + 2 V - 



dF a \ | d 2 F a \ df a 

dr\ ) drj 2 J drj 



{dP d^ F 1 

t4 + nl(6 + WaGr)^- + (2 + i Ma Gr)^| /„ 



dH 

-G a *—-?-sm6 + ifi a G a *H a cosO, 
drj 



(25) 



d 2 H n - dH r , 

= ~3Pr a F a - 



drf dr\ 



dF 1 

til + Pr a (-1 + i^ a Gr)^\ H a 



-Pr a /G a *(2 + in a Gr)-^f a , (26) 
where fi a = kSo is the dimensionless wave number normalized by the length So, and G a * = —dTa^/drj^—o, whose value 



depends on the Prandtl number. In the stability analysis, [15f v a — —dip a /dx and dT aif /dx were neglected because 
the derivatives of the unperturbed fields quantities F a and T a * with respect to x were assumed to be much smaller 
than those with respect to y. In this paper, however, these quantities in Eqs. (|25|) and (|26|) are retained because if 
we neglect them, and v p do not converge zero as [i a approaches zero. 



E. Equations of flow and temperature distributions in the water layer 

The stream function disturbance and temperature disturbance T[ in the water layer are assumed to be of the 
form: 0-[3 

ip'l = uiofi(y*)Ckexp[a-t + ikx], T( = i?;(y*)GCfcexp[ert + ikx], (27) 

where y* — y/ho, and /; and Hi are the dimensionless disturbance amplitude functions. It is also assumed that the 
unperturbed temperature distribution in the water layer is linear, then Gi = ^dTi/dy\ y= h — (T s i — Ti a )/ho. 
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When ipi — rpi is substituted into Eqs. ([T]) and ([2]), the perturbed part yields the following Orr-Sommerfeld 
equation for 0,0] 

= (2^ + inReVu) - [pf + mRei UlJu + ^) } /j, (28) 

where pi = kho is the dimensionless wave number normalized by the length ho, £//*(?/*) is the dimensionless velocity 
distribution in the water layer in the unperturbed state, and Rei = uioho/vi = 3Q/(2li>[) is the Reynolds number. 
Here Q/l is the water supply rate per width. 

When Ti —Ti + T[ are substituted into flSJ), the perturbed part yields the equation for Hi: 0,0] 

— - = itf + mPeiUi^Hi - i^Pet-^fi, (29) 

where Ti*(y*) = (T;(y*) — T s i)/(T s i — T/ Q ) = — is the dimensionless temperature distribution in the water layer in 
the unperturbed state, and Pe; = uiqIio/ki = 3Q/(21ki) is the Peclet number. 



F. Linearization of boundary conditions 



First, linearizing Eq. (jTU)) at y 



yields, to the first order in 
, dUu 



dfi_ 

dy* 



dy* 



0. 



/*L=o=0, 



(30) 



From the linearization of Eq. ((TTj) at y — ho, the relation between the amplitude of the water- air surface and that of 
the ice-water interface is obtained: = — (fi\ yt =i/U~i*\ yt =i)(k- 0-0] 
Second, linearizing Eq. (fl2|) at y — ho yields, to the zeroth order in 



d£a 
drj 



and to the first order in 



dfa 

drj 



d 2 F n 



T)=0 



drj 2 











F, 







U a 






uio 


{dU u 




(dh_ 


ho 


UaO 


\ dy* 


y,=i 


\dy* 



(31) 



fa 1 77=0 



fl\y„=l ) C^*|j/»=1 ( ■ 



UaO 



-Ul*\y, 



(32) 



The values of u a o are 0.38 m/s at x = 0.1 m and 1.2 m/s at x = 1.0 m for AT a ~ 10 °C. On the other hand, the 
surface velocity of the water layer, uio ~ [gs\n9/{2vi)] 1 ^[iQ/{2l)] 2 / 3 ', is about 0.78 ~ 3.62 cm/s for typical values 
of Q/l = 10 ~ 100 [(ml/h)/cm] and 9 = tt/2. It should be noted that the velocity of the water film flow is much 
less than that of airflow. Therefore, the first equation in (|3"TT) and the second equation in (|3"2")l are approximated as 
dFa/drjl^o = and / a |r/=o = ; respectively. Even though a thin fluid layer of water flows down the ice surface, the 
no-slip condition at the water-air surface of the flowing water film is nearly satisfied for the velocities: u a = dip a /dy, 
v a = —dtpa/dx and v' a = —dip' a /dx. The values of do are 3.7 mm at x = 0.1 m and 6.6 mm at x — 1.0 m for AT a = 10 
°C. On the other hand, the mean thickness of the water layer, 

ho = [2>vil{gsm6)Q/l] x /*, is about 53 - 115 /im for 
values of Q/l = 10 ~ 100 [(ml/h)/cm] and 9 = tt/2. Since So/ ho ^> 1, the second term on the right hand side of the 
first equation in Eq. (|3"2"|) cannot be neglected. Hence, the no-slip condition cannot be applied to u' a = dtp' a /dy at the 
water-air surface of the flowing water film. 

Third, linearizing Eq. (| 13[) at y = ho yields, to the zeroth order in 



dUu 



dy* 



PaV a (u a od 2 F a /dr] 2 \ n= o) /$o _ 



piviuio/ho 



— R r 



(33) 



and to the first order in 



d 2 f, 



dyl 

PaVa f h \ 2 U a ( d 2 f a 



d 2 U u 



pm \S J uio I drj 



V dy 2 
d 3 F a 



rj=o drf 



y,=n 

2 



?;=0 



UlJ\ yt= l + Pi J fl\y, = l 
+ Ptfa\ V =0 \ fl\y, = l/Uu\y, = l, 



(34) 
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where R Tal on the right hand side of Eq. (l33l) should be nearly considered as the ratio of the shear stress of 
airflow at the water-air surface to that of the water film flow at the ice-water interface. It is assumed that 
p[ = j o;wf o n;(y*)Cfe//ioexp[crt + ikx] and p' a — p a WQ n a (r/)^fc/^oexp[crt + ikx], where 11; and II a are dimensionless 
amplitudes. Substituting these forms into Eq. (|14[) and linearizing them at y = ho yields, to the first order in 



where 



d 3 fi 



{iHiRei(Ju\y m =i + 3/zf) j-^- 
cty* 



HiRei 



PaVg ( ho\ 3 UaO f d 3 fa 

Pin \S J u w \ dr] 3 
d 2 F, 



-ip a Gr 



drj 2 



r/=0 



77=0 

f a \n=o 4 



dy 

ijJLa.Gr 



dFa 

drj 



a/Uu\ y ,=i j fi\ y ,=i 

dfa 



77=0 



•3^ 



drj 



r/=0 



G a *H a \ v=0 sind > fi\y t —i/Ui»\ yil= i, 



(35) 



2(cat6)fit 



(36) 



represents a parameter relevant to the restoring force due to the surface tension and gravity acting on the water-air 
surface. 0, Q Here a = [7/ (/?/<?)] ^ 2 is the capillary length associated with the surface tension 7 of the water-air 
surface. (10| 

From the boundary condition (1331) and the no-slip condition Ui*\ y „=o = 0, the velocity profile in the water layer 
is given by i7;„ — yt + (R~r a i ~ 2)j/#. However, R Tal is extremely small because the ratio of the viscosity of air 
to that of water, p a v a /pivi, as well as Hq/Sq are much smaller than 1. Therefore, the shear stress-free condition, 
dUu / dy*\ V i,=i = 0, holds at the unperturbed water-air surface. Thus the velocity profile in the water layer is still the 
half-parabolic form, Uu — yl — 2y», so that the values of Ui*\ y ,=i = — 1, dUi*/ 'dy*\ Vr= o = —2 and d 2 Ui* / 'dy1\ Vr= i = 2 
are used in the above boundary conditions. Similarly, since p a v a i ' p\v\ <C 1 and ho/60 < 1 on the right hand side of 
Eqs. (jMJ) and (|33|) . the influence of the perturbed part of shear and normal stresses due to airflow on the water film 
flow at the water-air surface is negligible. Therefore, the boundary conditions for the shear and normal stresses at 
the perturbed water-air surface become the same as those used in the previous papers. 0-01 



Finally, linearizing Eq. (|17p at y — h yields, to the zeroth order in 
order in 



T, 



= -1 T 



Hl\y m =l + fl\y, = l/Uu\y t= l = 0, -ff a |?) = — L 

Linearizing Eq. (fTSj) at y = ho yields, to the first order in 



= 1, and to the first 



(37) 



dHj 




ho , 


' dH a 




dy* 


3/,=l 


do » 


\ dr] 


r7=0 ) 



fi\ y ,=i/Uu\ Vt= i — 0. 



It is convenient to define 




77=0 



Q'(i) 



ho / dm 



(i) 



So \ dr] 



7) = 



(38) 



(39) 



which represents the real and imaginary parts of the perturbed part of the air temperature gradient at the water-air 
surface. It should be noted that Eq. (|2"5)) can be independently solved with the boundary conditions (|3T)]) . and 
(|35f without considering the influence of airflow. Therefore, /; in Eqs. (l37l) and (|38|) is the same form as that in 
the absence of airflow. The perturbed part of temperature in the water layer is affected by the airflow through the 
perturbed part of the air temperature gradient in Eq. (|38|) . 



G. Dispersion relation 



From the perturbed part of Eqs. (jT5j) and (|16l) . the dispersion relation for the perturbation of the ice- water interface 
is given by 



v r m 

dy* 



+ Kf^ l (H l \ y , = o-l)\, 



(40) 



9 



where Kf = K s /Ki = 3.96 is the ratio of the thermal conductivity of ice to that of water. The real and imaginary 
parts of Eq. (j40|) give the dimensionless amplification rate = / (V / ho) and the dimensionless phase velocity 



v 



<r( % '/(kV), respectively, 



dH, 



(r) 



dy* 



+ Kfn(H<X, s0 -l), (41) 

3/„=0 



dH t 



(i) 

+ K!wHV\ ym=0 ] ) (42) 
y»=o I 



to \ dy 

where h\ and H^ are the real and imaginary parts of Hi. 

The numerical procedure for obtaining the wavelength and phase velocity of ice ripples is as follows. First, Eq. (|28|) 
is solved with the boundary conditions (jMf and (|3"5)) . Substituting the obtained solution f t into Eq. (|3"2"j) . then 
Eqs. (|2"2"j) . (j2"3")l . (|25l) and (j2"6")l must be solved simultaneously for a given Gr with the following boundary conditions: 
Eq. d3U), dFa/dT)\ n=00 = 0, T Q »|^ =0 = 1, T *|»j=oo = 0, Eq. (J35J), df a /dr/\ v=00 = 0, / |^ =00 = 0, H a \ n=0 = 1 and 
HJ v =oo = 0. Here, it is assumed that u' a \ y=oa = d%p' a / dy\ y=OQ = 0, v' a \ y=00 = —dty'J dx\ v=ao = 0, and T' a \ y=00 = 0. 
[15{ Substituting the obtained solutions /j and i? a into the boundary conditions ([3"Tf and (I3"5)) . Eq. (j2"9")l is solved. 
Finally, substituting the obtained solution Hi into Eqs. (|4ip and (|42p and replacing /i/ with (ho/8o)fj, a , it is possible 

(r) 

to calculate the amplification rate <r* and phase velocity u p * with respect to ^ a . 

III. RESULTS 

A. Solutions of temperature distributions in the air boundary layer 

In the absence of airflow, Eq. (|2"6")l yields d 2 H a /drj 2 = n 2 a H a . With the boundary conditions H a \ v —Q = 1 and 
-ffa|r)=oo = 0, the solution is given by H a — exp(— /U a ?y), and hence Ga^ — {ho/8o)fjt a — kha — and Ga^ — 0. In 
the presence of airflow, as shown in Fig. [2] (a), H^ decreases more rapidly than the exponential function, and Ha^ 
acquires non-zero values. Therefore, as shown in Fig. [5] (b), the value of Ga^ is greater than [i\ and Ga*' acquires 
non-zero values. 

From Eq. (|18[) . the energy conservation equation at the unperturbed water-air surface is — KidTi / dy\ y= h = 
—K a dT a /dy\ y= h - When the linear temperature profile Tj in the water layer and the exact temperature profile T a * 
in the air boundary layer are substituted into the above energy conservation equation, Ti a in Eq. (|17p is obtained as 

Tia ~ T s i + — ^ ° — T^. (43) 

Kl Oo/Ga* 

From Eq. (|16[) . the energy conservation equation at the unperturbed ice- water interface is LV = Ki(T s i — Ti a )/ho. 
Substituting Eq. ([4"3"f into this equation yields 

K T 

a _°° . (44) 

If <WG a * in Eqs. (j4*3"|) and (PHf is considered as <5 represented in Fig. [1] (a), then Ti a and V" in the previous papers 
[5H8J or V mentioned in the Introduction in this paper are obtained. The linear temperature profile in the air assumed 
in the previous papers is shown by the dashed line in Fig. [2] (a), which is expressed as T a * = 1 — G a *rj. Here, 
G a * = —dTa^/drjlri—o can be estimated numerically yielding a value of about 0.5. Using our notation, the boundary 
layer thickness in paper Q is expressed as S — C5o/V%, which must be equal to 6 = <$o/G a *. From this, the parameter 
G is determined as G = V2/G a * ~ 2.8. Since G a * is obtained from the solution of Eqs. (|22|) and ([23)) . S depends on 
the Prandtl number of air. 

B. Approximate solutions of flow and temperature distributions in the water layer 



Since 5q is of the same order as the characteristic length scale of ripples, we cannot use the long wavelength 
approximation, the higher order of fj, a in Eqs. ([23)1 and (f2"6"]l have to be retained. On the other hand, since the water 
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FIG. 2: For Q/l — 50 [(ml/h)/cm], 6 — it/2, x = 1.0 m and AT a = 10 °C, (a) air temperature distribution T a *, and distributions 
of exp(— HaTj), Ha and Ha^ at the dimensionless wave number of fx a = 4.8. (b) perturbed part of air temperature gradient 
G' a = ho/6o(— dHa/dr/lrf^o) at the water-air surface: in the absence of airflow G' a — fxi; in the presence of airflow G' a and Gl' 1 ' 
are the real and imaginary parts of G' a . Here /j a = 10 corresponds to the wavelength of 4.1 mm when So = 6.6 mm. 



layer thickness ho is much less than the characteristic length scale of ripples, we can neglect the higher order of /i; in 
Eqs. l[28|). (|29| . (|34| and ([35]). Using the long wavelength approximation, /; and Hi can be calculated approximately 
as in the previous papers. [J-Q 

Transferring the variable to z = 1 — j/*, the general solution of (f29")l is expressed as: 0,0] 

H l {z)=C l( j> l {z) + C 2 ^{z)+mPei f {(f> 2 {z)Mz') ~ Mz)Mz')} fi(z')dz' , (45) 



where 0i and </>2 are solutions of the homogeneous equation From Eqs. (j3"7| and (j3"8")l . we obtain Ci = /i| z =o and 
C 2 = /ioAo(-d-ffa/dr?|,=o)/jU=o, respectively, because <pi\ z=0 = 1, 02 U=o = 0, d0i/dz| z=o = and d(j> 2 /dz\ z=Q = 1. 
Consequently, 77; is expressed as 



So V rfr 7 



?;=0 



+iViPei {fa^M*') - faWfaiz')} fi{z')dz' '. (46) 



For typical values of ho and uio, Rei ~ 1 and Pe; ~ 10; then fiiRei -C 1 and /J-iPei ~ 1 for the length scale of ripples 
on icicles. Therefore, we can neglect the niRei term in Eqs. ([28]) and (|35|) . This corresponds to neglecting the inertia 
term of the full Orr-Sommerfeld equation. Q Furthermore, the expansion of 4>i and 4> 2 with respect to [iiP&i up to 
the first order is sufficient. Indeed, the justification for these approximations was confirmed by our recent numerical 
analysis. Q Hence, it is sufficient to use the following approximate solutions: 

fl(z) = —}—(6 + iaz-6z 2 -iaz 3 ), (47) 
b + ta 

Mz) = l-i( K \z 2 --^z^iPei, (48) 

Mz)=z-iQ:z 3 -^z^ l Pe l . (49) 

Since the direction of the x axis in Fig. [T]is opposite to that in the previous papers, [343 we n °te that the sign of 
in this paper is opposite. This leads to different functional forms of <j>\ and <f> 2 from those in the previous papers. 
In the presence of airflow, using the approximate solutions (j47|) . (|48| and (J49]) . Eqs. (jS} and (|42|) yield 

(r) _ Ga (r) {36 - f qfciPeQ} + G'^ {6a + 9^Pe;} - f afctPeQ 
ff * ~ 36 + c? 

_ G'P {36 - ^a( W Pe,)} + G£ w {6a + - ^(w^) - a 2 

+*i « , (50) 
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1 



\a 2 {^P ei ) + Ga [r) {6a + 9 W Pe,} - gT ] {36 - §a( W Pe;)} 



/(i) 



36 + a 2 

6a - X a 2( WjPe;) + Q'iT) | 6a + » w p e ,} - G' Q W {36 - ^a( M; Pe,)} 



36 + a 2 



(51) 



On the other hand, in the absence of airflow, since G'^ = \i\ and G'^ = as mentioned above, Eqs. (|50l) and (f5Tj) 
reduce to the previous dispersion relation: [j, 0] 



(r) _ m {36 - |a(^/Pe ; )} - ^a(piPei) fit {36 - ^a(mPei)} - ^a(mPei) - a 2 
C* — — - h -fV; /i; 



36 + a 2 



36 + a 2 



(52) 



Mi 



-\a 2 (mPei)+ fii{6a + 9mPei} _ 6a - ^a 2 (/i ; Pe ; ) + w {6a + f mPei) 

h A ; ^; 



36 + a 2 



36 + a 2 



(53) 



C. Wavelength and translation velocity of ripples 



For the water supply rate per width Q/l = 50 [(ml/h)/cm] and the angle 9 = tt/2, Figs. [3] (a) and (b) show 
numerically obtained the dimensionless amplification rate o* = a^/(V/ho) and the dimensionless translation 
velocity v plf = v p /V versus dimensionless wave number fi a = kSo, respectively. The wave number of ripples that one 
expects to observe is that for which the amplification rate is the maximum. We also define the value of w p * from the 

wave number at which er* acquires a maximum value. In the presence of airflow, a), acquires a maximum value of 

a (r) 
u *max 



r M = 

= 0.59 as represented by the dashed line in Fig. 



acquires a maximum value of crimax = 0.054 at /i a = 4.3 (dashed line in Fig. [3] (a)), 



0.085 at [i a = 4.8 (solid line in Fig. [3] (a)). Since the wave number fc is normalized by Sq, the corresponding 
wavelength is 8.6 mm from A = 2ttSq/ fi a - Here we have used So = 6.6 mm estimated from the two parameters x = 1.0 
m and AT a = 10 °C. At ix a = 4.8, v ple = 0.48 as represented by the solid line in Fig. [3] (b). On the other hand, in 
the absence of airflow, Eq. ([52 

which corresponds to the wavelength of A = 9.6mm. At \i a = 4.3 
E](b). 

Any disturbance near the solidification front can be initiated by non-uniformity in temperature in the vicinity of 
the ice-water interface. Since the water layer considered here is very thin, we cannot neglect the influence of external 
disturbance at the water-air surface on the growth condition of the ice-water interface. In order to determine the 
growth condition from the dispersion relation (|40p . it is necessary to obtain the perturbed temperature amplitude 
Hi in the water layer. Hi must satisfy the boundary condition ([3"5]) which includes the perturbed air temperature 
gradient at the water-air surface. Using Eq. ([39)1 and Ui*\ y ,=i — —1, the real and imaginary parts of Eq. ([38)) can be 
written as follows: 



dH 



(r) 



dy* 



Ga r) fi 



Mi 



GTfi 



Since Ga = ^i and G a — in the absence of airflow, Eq. 



dH 



(i) 



dy* 



G , W/ ( 



Mi 



(54) 



54]) reduces to the previous results: Q 



dH, 



M 



#Mi 



dH) 



dy* 



(55) 



The solid and dashed lines in Fig. [3] (c) show the behaviour of —dH^ 



«// P> l« 



_! in Eq. ([54]) and of 

and vifi r) \y, 



f ) /dy m \ v . =1 ,-dH®/dy m \ v .-_ 

u A*/// lj/»=i m Eq. ([55[) with respect to fi a . It can be seen that — dH^ J dy*\ y ^—\ ana /j; /, , 7: . i 
for small /i a . In the absence of airflow, the rate of latent heat loss due to thermal diffusion from the water- air surface 
to the air changes locally by the water-air surface disturbance. On the other hand, in the presence of airflow, the 
rate of latent heat loss is enhanced by the airflow, more so than in the case of thermal diffusion. However, as shown 
in Fig. [3] (c), non- uniformity of the rate of latent heat loss at the water-air surface decreases with an increase in 
Ha because of the action of the restoring force on the water-air surface, which causes the amplitude of the water-air 
surface disturbance to decrease. 0-0 This effect is due to the parameter a in /; in Eqs. ([54")) and (|55[) and is more 
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FIG. 3: For Q/l = 50 [(ml/h)/cm], 8 — ty/2 and 8q = 6.6 mm, (a) dimensionless amplification rate ai = cr' r ' /(V /ho) versus 
dimensionless wave number jj, a = kSo; (b) dimensionless phase velocity v pt — v p /V versus dimensionless wave number fi a - Solid 
and dashed lines indicate the presence and absence of airflow, respectively, (c) The behaviour of the real and imaginary parts 
of the perturbed temperature gradient at the water-air surface with respect to \x a , in the presence of airflow (solid lines) and 
in the absence of airflow (dashed lines). Here \x a = 100 corresponds to the wavelength of 413 fim when So = 6.6 mm. 



effective for large wave numbers. The physical meaning that the values of —dH^/dy^y^^ in Eq. (f54"f and Hif^ |y„=i 
in Eq. (|55|) are not zero will be discussed in IIIIDI 

An approximation of Eq. (|50|) makes the above discussion more clear. We note that the second term in Eq. (|50|) 

(r) 

is smaller than the first term, and the wave number at which cr* acquires a maximum value is almost the same as 
that without the second term. Q Therefore, extracting the most dominant term from the first term in (|50p and using 
we obtain 

Jr) _ 36Gq (r) - ZgjmPet) _ ^ f(r) Pe t f a\ 2 4 

36^ " Gq ~T2~UJ ^' ( j 

at 8 = tt/2. As mentioned above, the non- uniformity of the air temperature gradient at the water-air surface is the 
trigger of the ice- water interface instability, which is represented by the positive term Ga r ' in Eq. (|56p . In the absence 
of airflow, since Ga^ = fj,i, we find from da^ /dfii = that <r^ acquires a maximum value at fii — [3(ho/ a) 2 / Pei] 1 / 3 . 
From this, an approximate formula is obtained to determine the wavelength of the ripples: A = 2ir(a 2 hoPei/2>) 1 / 3 , 
16, 7\ as mentioned in the Introduction in this paper. On the other hand, in the presence of airflow, the value of 

Ga^ is greater than /i;, as shown in Fig. [2] (b). This indicates that the natural convection airflow enhances the 
destabilization of the ice-water interface compared to the destabilization due to the thermal diffusion. However, it is 

l(r) 

difficult to express the dependence of G a on /i a analytically. The stabilization of the ice- water interface is dominated 
by the negative term in Eq. (|56[) . The stabilization mechanism due to the action of the restoring force of the surface 

tension and gravity on the water-air surface is not relevant to the airflow. Although the value of c^max m the presence 
of airflow is greater than that in its absence, the wavelengths determined from the most unstable mode have nearly 
the same value in both cases. However, there is a considerable difference in v plf . In the absence of airflow, u p * > 



13 



for all fi a , as shown by the dashed line in Fig. [3] (b). On the other hand, in the presence of airflow, u p * has negative 

values for a small wave number region because the terms with Ga^ in Eq. (|5Tj) are the most dominant. The solid line 
in Fig. [3] (b) indicates that the sign of u p * changes from negative to positive at fi a = 3.7. What determines the sign 
of v p * will be discussed in lHIDI 




0.2 i \ ; ; - 

o L i i i i 

20 40 60 80 100 

(c) 0/1 [ (m I /h) /cm] 

FIG. 4: (a) The wavelength versus sin# at Q/l =160/3 [(ml/h)/cm]. (b) The wavelength versus Q/l at 9 = n/2. (c) The phase 
velocity versus Q/l at 8 — tt/2. Solid lines indicate the absence of airflow; [3,0] dashed lines and dashed-dotted lines indicate 
the presence of airflow for Gr = 609 and Gr = 108. 

Figures E] (a) and (b) show the dependence of the wavelength of ripples on sin# at Q/l — 160/3 [(ml/h)/cm] and 
that on Q/l [(ml/h)/cm] at 9 = it/ 2, respectively. We have determined these wavelengths from the value of \i a at 

(r) 

which <ri acquires a maximum value for a given Q/l and 6. The solid lines indicate the absence of airflow, whereas 
the dashed and dashed-dotted lines indicate the presence of airflow. Although the wavelengths of ripples in the 
presence of airflow are slightly shorter than those in the absence of airflow, the dependence of the wavelengths on the 
angles and water supply rates shows almost the same behaviour as the experimental results (see Figure 8 Q)- When 
determining the wavelengths in Fig. |4](a) and (b), we have used 5q = 6.6 mm (x = 1.0 m and AT a = 10 °C) for the 
dashed lines and 5o = 3.7 mm (x = 0.1 m and AT a = 10 °C) for the dashed-dotted lines. 

Table H] shows the wavelengths obtained from various values of Sq — Ax/Gr using different combination of x and 
AT a . It is found that the wavelength increases with the increase of both x and AT a . This suggests that G'^ in Eq. 
(l56l) must include the modified local Grashof number Gr. However, the dependence of the wavelength A on i and 
AT a in Gr is extremely small compared to that of V, Ti a and S. This result is relevant to the fact that the wavelength 
of ripples on icicles is nearly independent of the vertical position of icicles and ambient air temperature. Table H] also 
shows that the value of v pse increases with a decrease in AT a and a decrease in x. Since v p * has positive values, the 
ripple with the most unstable mode moves only upwards. Figure [4] (c) shows the dependence of v p * on Q/l. The 
range of variation of v pif on Q/l in the presence of airflow (dashed and dashed-dotted lines) is larger than that in 
the absence of airflow (solid line), and the dependence of v pit on Gr is larger than that of the wavelength A on Gr. 
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Therefore, we can say that i> p * is sensitive to the parameters characterizing the air boundary layer. 



TABLE I: For Q/l = 50 [(ml/h)/cm] and 8 — it/2, the dependence of modified local Grashof number, Gr, ice growth rate, 
V, temperature at water-air surface, T; a , thickness of thermal boundary layer, 5 = 5o/G a *, wavelength of ripple, A, and 
dimensionless translation velocity of ripple, v p *, on air temperature far away, Too, and position from the bottom of the gutter, 
x. 



x = 1.0 m 



Too (°C) 


Gr 


V (mm/h) 


Ti a (xlO- 3 °C) 


8 (mm) 


A (mm) 


v p * 


-5 


512 


0.08 


-1.2 


16.1 


8.3 


0.65 


-10 


609 


0.20 


-2.9 


13.4 


8.6 


0.48 


-15 


674 


0.32 


-4.9 


12.1 


8.7 


0.41 


-20 


724 


0.47 


-7.0 


11.2 


8.7 


0.37 


x = 0.1 m 


-5 


91 


0.14 


-2.1 


9.6 


7.9 


0.89 


-10 


108 


0.33 


-5.0 


7.9 


8.3 


0.71 


-15 


120 


0.56 


-8.4 


7.0 


8.1 


0.63 


-20 


129 


0.81 


-12.1 


6.5 


8.5 


0.56 



D. Heat flux at the ice- water interface and water-air surface 



We assume a dimensionless small perturbation of the ice-water interface with an infinitesimal initial amplitude 
5b = (k/h : 



y* — C* = <J&Im[exp(<7*i» + i/^x*)] = 8b(t*) sin[/z/(a;* — u p *t*)], 



(57) 



where er* = a/(V /ho), t* = (V/ho)t, x* = x/ho, Sb{t*) = #6exp(<7* i*) and Im denotes the imaginary part of its 
argument. The corresponding perturbation of the water-air surface with an infinitesimal initial amplitude St — £k/ho 
is given by 



y* = £* = 1 + Im[<5 t exp(o\T* + ifux*)] 

= 1 + [(// r) |y,=i) 2 + (fi l) \y,=i) 2 } 1/2 Sb(t.) sin[ W (z* - v p M) - 6 5 J, 



(58) 



where the relation 5t — fi\ y ,=iSb for the amplitude is used, and 0£„ is a phase difference between the water-air surface 
and the ice-water interface. Since //| y »=i depends on the wave number through the parameter a, the amplitude and 
phase of the water-air surface relative to the ice- water interface change depending on the wavelength of the ice- water 
interface disturbance. |4|-|7| 

The temperatures in the water layer and ice are expressed in the dimensionless forms: 0] 



Ti*{y*) = Tl ^*^ r^ sl = -y* + 5 6 Im[i?i(2/*)exp(<T»t* + i/ija;*)], 

± si — ± la 



(59) 



T s *(y*) = Ts ^*^ = 5bexp(jj,iy*)Im[(Ht\y m =o - l)exp(cr*t* + ifnx*)] 



T s l — Ti a 

and the temperature in the air boundary layer is expressed as 

df, 



T a *{rj) = T a *(rj) +Im 



dij 



v=oj o 



(60) 



(61) 



where we note that y is normalized by ho in the water layer and ice, but y is normalized by $o m the air boundary 
layer. 

Wc define the perturbed part of dimensionless heat flux from the ice-water interface to the water and from the ice 
to the ice-water interface, as qi* = Im[— 9T i ' <t /9y»| 2/ , = ^,] and q s * = Im[— KfdT^/dy*\ yr= ^], respectively, where 
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and T' s1f represent the perturbed terms in Eqs. 
to the water and ice is expressed as follows: |7j 



59l and (l60l). Hence, the total heat flux from the ice- water interface 



qi* 



<5jJm 



dHi 

dy* 



dH, 



(r) 



y«=o 



xf w (ff/ r) L=o-i) 



fff W (-ffily.=o 

\ 2 / 



1) > exp(cr*i* + imx*) 



y*=o 



dH, 



(i) 



dy* 



1/2 



Xi5 6 (t*) sin[^(:c* - - O g!s „], 



(62) 



where 9l3t is a phase difference between the total heat flux qi s * at y* = £* and the ice-water interface. We also define 
the perturbed part of dimensionless heat flux from the water-air surface to the air as q a * = lm[—dT ar /drj\ v= ^/s a ], 
where T' a ^ = T' a /(Ti a — T^) represents the perturbed term in Eq. (jST]) . Hence, 



Cfc T 

-— 1m 



dH n 



))=0 



/;ls/„=i ex P( ( M* + ipi x *) 



dij 

(<# p) /, W lv.=i 
x £);,(£*) sin[/i;(a;* — w P *i*) — © 9q , 



r"(0 /W| 



1/2 



(63) 



where Oq a , is a phase difference between the heat flux q a * at rj — £' /So and the ice-water interface. 

Figures [5] (a) and (b) illustrate the time evolution of the ice- water interface with an initial amplitude of 8b — 0.05, 
for the wave number /x a = 4.3 in the absence of airflow and for [i a = 4.8 in the presence of airflow, respectively. The 

(r) 

respective wave number represents the fastest growing mode, at which ct; acquires a maximum value, as shown by the 
dashed and solid lines in Fig. [3] (a) . The arrows on the ice- water interface and the water-air surface show the position 
of the maximum of qi s * and that of q aSr 
(-dH^/dy4 y , =1 ) 2 }^ 2 S b (t,)sm[^{x, - 



Using Eq. ([54| . Eq. (|63|) can be written as q a * = [{—dH( r '/dy, 



— 9a J. Therefore, non-zero values of 



-dH^/dy* 



,=l in Eq. (|541 



contribute to the imaginary part of q a *, and cause the phase shift of q a * relative to the ice- water interface. 

In the absence of airflow, as shown in Fig. [S](a), q a * is largest at each protruded part of the water-air surface because 
the isotherm in the air is symmetrical around the protruded part. As shown in Fig. [5] (c), the water-air surface shifts 

to the positive X* direction by 0£„ relative to the ice-water interface. In the absence of airflow, since Ga 



'(r) 



and 



Ga 



'(<) 



0, Eq. 



yields q a * = \y,=i) 2 + (/; \ y ,=i) 2 ] 1/2 5b{t*) smlfiitx* - w p *i*) - 6 9a J. Comparing this to 



Eq. (f58|) . it is found that the position of the maximum of q a * also shifts to the positive direction by 9ai = 6^,. 
However, the position of the maximum of qi s * shifts by gis „ to the upper side of the protruded part of the ice-water 
interface. 

On the other hand, in the presence of an upward airflow shown by the dashed arrow in Fig. [5] (b), the isotherms 
in the air boundary layer are no longer symmetrical around each protruded part. The isotherms become closer on 
the lower side of the protruded part of the water-air surface due to the upward airflow. Hence, q a * is largest on the 
lower side of the protruded part, as shown in Fig. [S] (b). By comparing Fig. [S] (a) to (b), first, it is found that the 
position of the maximum of q a * in the absence of airflow is always on the protruded part of the water-air surface, but 
that position changes by the presence of airflow and depends on the wave number fi a . As shown in Fig. [S] (c), the 
sign of Oq at in the presence of airflow changes from negative to positive value at /i a ~ 5.7, which corresponds to the 

change of sign of —dH[ l ^/dy^\y t= i in Fig. |3](c). Second, there is a critical difference between the phase shift Q qist in 
the absence of airflow and that in its presence. In the absence of airflow, the position of the maximum of qi s * shifts 
to the upper side of the protruded part of the ice-water interface with an increase in fj, a (see Q qist (no airflow) in Fig. 
[5] (c)). In this case, the sign of u p * is positive as shown by the dashed line in Fig. [3] (b). Figure [S] (a) shows that 
the ripple at fi a = 4.3 moves upwards at w p * = 0.59. The displacement in the dimensional form is about 11 ho after 

the dimensionless time l/c*max = 1/0.054. On the other hand, in the presence of upward airflow, the position of the 
maximum of qi s „ is on the lower side of the protruded part of the ice-water interface for < fi a < 3.7, whereas that 
is on the upper side for [i a > 3.7 (see Q qisr (airflow) in Fig. [5] (c)). We showed that the sign of changes from 
negative to positive at fi a = 3.7 by the solid line in Fig. [3](b). Therefore, the sign of v p * is related to the sign of Q qis , ■ 
The ripples move down in the mode < fi a < 3.7, whereas they move up in the mode /x a > 3.7. However, the ripple 
with the most unstable mode of /i a = 4.8 is expected to be observed. Figure [5] (b) shows that the ripple at [i a = 4.8 
moves upwards at v p * — 0.48. The displacement in the dimensional form is about 6 ho after the dimensionless time 

l/o&L = 1/0.085. 
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FIG. 5: (a) and (b) are illustrations of the time evolution of an initial disturbance of the ice-water interface from t* = to 
U = The solid arrows in the water film and the dashed arrows in the air boundary layer show the direction of the 

supercooled water flow and airflow, respectively. The arrows attached qi a * and q a * are the maximum point of heat flux at the 
ice-water interface and water-air surface, respectively, (a) represents the disturbance of /j, a = 4.3 in the absence of airflow, (b) 
represents the disturbance of fi a — 4.8 in the presence of airflow, (c) represents the phase shift of the water-air surface, Q^ kt , 
total heat flux at the ice-water interface, & qist , and heat flux at the water-air surface, & Qa , , relative to the ice-water interface 
with respect to fj, a . 



Substituting Eqs. ([55]) and (JBOJ) into Eq. ([L3|). the dimensionless form of AT s i can be written as: 

AT sl * = [{H[ r) \ y , =0 - l) 2 + (i? ; W | y , =0 ) 2 ] 1/2 (5fc(t*)sin[ A1/ (a; >t - v p M) - 9 T J, (64) 

where 0*r Cl is a phase difference between the temperature at = £* and the ice-water interface. H{ and in (|64|) 
are determined by solving with the boundary conditions ([57)1 and . Since the water film flow is not affected by 
the natural convection airflow, the forms of Uu\ y „ = \ and fi\ y , = i in (|37p and are the same as those in the absence 
of airflow. However, dH a /dj^\ v= o in (|38[) in the presence of airflow is different from that in the absence of airflow. 
As a result of the change in the temperature gradient at the water-air surface, the solution Hi changes and causes 
different distribution of AT S ;, at the ice-water interface. The position of the maximum of AT S /* changes depending 
on that of the rate of latent heat loss at the water-air surface. 0] That is why AT S ; in Eq. (fT5)) was considered as 
the spatial temperature non-uniformity caused by the external disturbance at the water-air surface. Heat flux q s * in 
the ice in the vicinity of the ice- water interface is caused due to the deviation AT S ;*, which contributes to the second 
term in Eq. (gOJ) • 
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IV. SUMMARY AND DISCUSSION 

A morphological instability theory has been elaborated for ice growth under a water film flow with a free surface 
and a natural convection airflow, within a linear stability analysis. This theory proposes a synthetic treatment of 
the heat flow in ice, water and air through a disturbed ice-water interface and water-air surface, thin water film 
flow and airflow, taking into account the influence of the shape of the water-air surface on the growth condition of 
the ice-water interface disturbance. Even though the natural convection airflow was introduced, the shear stress-free 
condition at the unperturbed water-air surface still held. Moreover, the influence of the perturbed part of shear and 
normal stresses due to natural convection airflow on the water film flow was negligible. Consequently, the perturbed 
distribution of water film flow could be obtained without considering the influence of the airflow. However, since the 
rate of latent heat loss from the water-air surface to the surrounding air is affected by the airflow, the perturbed 
temperature distribution in the water layer is different from that in the absence of airflow. In the absence of airflow, 
the position of the maximum of heat flux q ast at the water-air surface is at the protruded part of the water-air surface. 
In the presence of airflow, that of q a * is not necessarily at the protruded part. Depending on the position of the 
maximum of g a *, that of <7/ s * at the ice- water interface changes. We find that the position of the maximum growth 
rate of the ice- water interface disturbance is shifted upward relative to the position of the maximum of q a * . We also 
find that although the airflow causes the amplification rate of the ice-water interface disturbance to increase by the 
enhancement of the rate of latent heat loss from the water-air surface to the surrounding air, the wavelength of ice 
ripples is not significantly affected by the natural convection airflow. On the other hand, the mean ice growth rate V 
and the ripple translation velocity v p depend on the parameters characterizing the air boundary layer. 

We mention the importance of the influence of the temperature distribution in water film flow on the growth 
condition of the ice-water interface disturbance even though the water layer is very thin. If we can neglect 
the temperature distribution within the water layer, and focus on only the temperature distribution in the air, 
Eq. (Tip]) is replaced by L{V + d(/dt) = —K a dT a /dy\ y= ^. Linearizing this equation at y = ho yields, to 
the zeroth order in V = — KaTcx/ (L8o/G a *), which is identical to Eq. (|4"4")) . The first order in gives 
a = (V/ho)(ho/So)(—dH a /dr]\ ri= o)fi\y :t= i, whose real part is approximately expressed as a* — Ga fi \y,=i — 
G'^ } ti l) \y,=i = (36Ga (r) + 6aG' Q W )/(36 + a 2 ) « G' Q (r) . Comparing this to Eq. (J56]), it is found from Fig. [2] (b) 
that the ice- water interface is always unstable because the stabilizing term is absent. It should be noted that the 
stabilizing term in Eq. (1551) was obtained from the solution of the perturbed temperature distribution in the water 
film flow. Although the heat transfer through the air boundary layer is the deciding factor in the growth rate V, in 
order to obtain the growth condition of the ice- water interface disturbance, it is important to determine the perturbed 
temperature distribution in the water layer as well as that in the air boundary layer. 

Although the wavelengths theoretically obtained in Table Q] are in agreement with experimental results, 0, [H, 0] 
several questions arise for the values of V and T/ a . The measured growth rates of icicle radius in the experiment was 
about 1.4 ~ 5.3 mm/h for the air temperatures range of —4.9 ~ —28.8 °C in the case of the zero wind speed. [l[ 
Also, the mean radial growth rate of ice grown on a 6-mm diameter round stick was 1.7 mm/h (see Figure 9 (a) 0). 
This experiment was conducted in a cold room, where large temperature fluctuations of ±3 °C around —9 °C were 
observed. Substituting the measured value into the energy conservation equation LV — — K{Fi a / 'ho at the ice-water 
interface, the degree of supercooling of the water layer is T[ a = —0.03 °C. Certainly, the values of T\ a and V calculated 
from Eqs. ([4*5]) and (f4"4l are less than the measured experimental values by one order of magnitude. If the value of 
the boundary layer thickness 8 is less than that estimated from the natural convection boundary layer, Eqs. (|43[) and 
(|44| suggest that the values of Ti a and V should increase. It is known that it is somewhat difficult to grow icicles 
with significant ripples in the steady calm conditions of icicle formation. [l[ Therefore, instead of assuming a calm 
environment for ice growth, different heat transfer mechanisms needs to be considered. Also, it is necessary to predict 
or to measure the mean ice growth rate V accurately in order to estimate the displacement of ripples. We have to 
be careful in the measurement of the displacement of ripples because v p depends on environmental conditions, as 
mentioned above. 

Finally, limitations of the proposed theory must be mentioned. First, it was assumed that ice was grown in a 
flat gutter on an inclined plane, considering a perturbation around the flat ice surface. However, as shown in Table 
HI for a given air temperature X^, since the ice growth rate V depends on x, the actual grown ice thickness on 
the gutter varies locally. If heat conduction through the ice to the substrate is negligible, the ice thickness bo in 
the unperturbed state is proportional to the time. Q The angle that tangent vector to the ice-water interface at 
(x,bo) makes with respect to the positive x direction is given by (f>(x,t) = cos _1 [{l + (dbo/dx) 2 }^ 1 ^ 2 ]. Making 
use of Eq. (|44|) . x = hox* and t = (ho/V)t*, (f>(x,t) gradually changes in time from the initial flat ice surface by 
cos _1 [{l + {(dV /dx)t} 2 } -1 / 2 ] = cos _1 [{l + {t»/(4a;»)} 2 } _1 / 2 ]. However, the change is negligible except for small 
because <C 1 even after 10 hours in the range of 0.1 < x < 1 m. The actual geometry of the icicle is that of 

an elongated carrot shape. Q In this case too, since icicle's surfaces are nearly vertical, we can neglect the change in 
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the slope dbo/dx in <f)(x,i) except for the tip region. Hence, the use of air boundary layer under the assumption of a 
flat ice surface is valid, [8J and the local variation in the thickness h and the surface velocity u; of the water film 
in the unperturbed state is negligible as ice grows. However, for ice growth on aircraft wings and aerial cables, the 
local change in (f> in time is remarkable compared to the icicle growth, so that we have to consider a morphological 
instability around curved ice surfaces in the unperturbed state, and ho and uiq are no longer constant over the curved 
ice surface. This is relevant to the problems on solidification on surfaces of arbitrary curvature. [l7[ 

Second, in our linear stability analysis, a small perturbation of the ice-water interface was assumed: y» = C* — 
5b (i*) sin[/Lti(x* — t> p *i*)]. However, since the amplitude <5&(t*) = <5i,exp(cr* <*) in £* and in the corresponding fields 

(r) 

increases exponentially with time when <r* > 0, the non-linear terms for the perturbation in the governing equations 
and boundary conditions are no longer small. Even though the linear approximation only describes the initial evolution 
of the perturbation, there was good agreement between the wavelengths predicted from our linear stability analysis 
and experimentally observed wavelengths of finite amplitude ripples. However, it is needless to say that the linear 
theory is unable to clarify further features related to ripple development, and the question arises of the value of the 
saturation amplitude, and of how the perturbation amplitude evolves towards this value. [l4| This leads us to extend 
the linear perturbation calculation to higher orders in the perturbation amplitude. [l8| Such an amplitude expansion 

generalizes the time evolution equation of the amplitude of the ice- water interface from d8b{t*) / 'dt„ — ai r ^Sb(t r ) to a 
nonlinear amplitude evolution equation. In order to implement it, algebraically complicated calculations are needed. 

Third, for the relatively weak flow considered here, the free shear stress condition at the water-air surface was still 
satisfied, and water film flow was driven by gravity only. However, in the presence of a strong airflow around aircraft 
wings and aerial cables, the water film flow is driven by gravity and aerodynamic forces. Due to strong air shear stress 
exerted on the water-sir surface, the distribution of water film flow must be modified from the half-parabolic form 
Ui* = yl — 2y* to Uu = vl + {Rr a i — 2)y», as discussed in III Fl It is known that the aerodynamic forces, as modified 
by the accreted ice, are significant in determining the wind drag and lift on iced structures. However, the traditional 
approach in wet icing modelings has been based on the mass and energy conservations only and have ignored the 
dynamics of the surface flow of unfrozen water. [l9| When airflow and water film flow are coupled, the distribution of 
shear and normal stresses at the water-air surface may influence the temperature distribution in the water layer. The 
action of an aerodynamic force on the water-air surface, and the resulting morphological instability of the ice-water 
interface have to be considered. These issues are beyond the scope of the analysis developed here. Removing these 
restrictions will be the subject of future research. 
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